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We use extensive classical molecular dynamics simulations 
to calculate the thermal conductivity of a model silica glass. 
Apart from the potential parameters, this is done with no 
other adjustable quantity and the standard equations of heat 
transport are used directly in the simulation box. The calcu- 
lations have been done between 10 and 1000 Kelvin and the 
results are in good agreement with the experimental data at 
temperatures above 20K. The plateau observed around lOK 
can be accounted for by correcting our results taking into 
account finite size effects in a phenomenological way. 

PACS numbers: 61.43.Fs, 61.20.Ja, 66.70.+f, 65.40. +g 



I. INTRODUCTION 

The thermal properties of glasses exhibit some specific 
and unusual features which are well known for quite some 
time ^ . These features are apparent in the specific heat 
and the thermal conductivity but we would like to focus 
here on the thermal conductivity k. The temperature 
dependence of k(T) can be separated in 3 distinct tem- 
perature domains: 

- At very low temperature {T < IK) the thermal conduc- 
tivity increases like T^. This increase can be explained 
within the tunneling model Q which has been proposed 
almost thirty years ago. 

- At intermediate temperatures (2 < T < 20K) the ther- 
mal conductivity exhibits a "plateau" for which several 
explanations have been given An extension of the 
tunneling model, the soft-potential model, has been pro- 
posed and gives a coherent description of the plateau by 
introducing the concept of "soft vibrations" 

- At high temperature, (T > 30K), k{T) rises smoothly 
and seems to saturate to a limiting value Kqo unlike crys- 
tals where k{T) ~ 1/T at elevated temperature. Re- 
cently this second rise of the thermal conductivity has 
also been explained within the soft-potential model ||] 
which appears to be able to account for all the thermal 
anomalies of glasses over the whole temperature range. 
Our aim here is not to propose a new or alternative expla- 
nation of the above mentioned anomalies. The purpose 
is to perform a molecular dynamics (MD) simulation on 
a model silica glass using a very widely used interaction 
potential (the so-called "BKS" potential [Q) without any 
pre-conception of the model able to explain the thermal 
anomalies of silica. This means that we do not add or 



inject an a priori quantity in the potential to reproduce 
a specific model. We use the standard definition of the 
heat transport coefficients that we calculate directly in 
our simulation box. In fact we introduce artificially in- 
side the system a "hot" and a "cold" plate which there- 
fore induce a heat flux. This flux creates a temperature 
gradient and once the steady state has been reached we 
can determine the thermal conductivity. By using plates 
compatible with the periodic boundary conditions we are 
able to calculate the thermal conductivity directly dur- 
ing the simulations without any additional parameter. 
This technique has been inspired by earlier studies |^ 
in which the plates were treated like hard walls and has 
mainly been applied to the calculation of the thermal 
conductivity in 1- or 2-dimensional systems |9pC[|. Nev- 
ertheless very recently Oligschleger and Schon applied 
the same method in a study of heat transport phenom- 
ena in crystalline and glassy samples (mainly selenium) 
In parallel to these studies which can be called in 
situ, other methods relying on the use of the density and 
heat flux correlation functio ns [[l2[ or on the Kubo and 
Greenwood-Kubo formalism |13| have been developed in 
order to determine the thermal conductivity of solids. 
Our results for the thermal conductivity obtained with 
the BKS potential compare reasonably well with the ex- 
perimental data. First of all the order of magnitude is 
correct above 20K and, at least in the range 20K-400K, 
a nice quantitative agreement is obtained. Furthermore, 
by taking care of finite-size corrections in a very sim- 
ple phenomenological way, we are able to reproduce the 
plateau around lOK. Of course, the very low tempera- 
ture behavior, which is known to be due to quantum 
effects, is out of the scope of such a classical calculation. 
This paper is organized in the following way. In section 
II we describe the modus operandi we have used to ob- 
tain the thermal conductivity. In section III we present 
first the results obtained directly from the MD simula- 
tions. Then we show the effect of finite-size corrections 
on these results and discuss our findings. In section IV 
we draw the major conclusions. 



II. MODUS OPERANDI 

Except the determination of k{T), the simulations are 
standard classical MD calculations on a microcanonical 
ensemble of 648 particles (216 Si02 molecules) interact- 
ing via the BKS potential. Like in a previous study 
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Q the particles are packed in a cubic box of edge 
length L = 21. 48 A (the density is approximately equal 
to 2.18g/cm'^) on which periodic boundary conditions are 
applied to simulate a macroscopic sample. The equations 
of motion are integrated using a fourth order Runge- 
Kutta algorithm with a time step At equal to 0.7fs. The 
glassy samples are obtained after a quench from the liq- 
uid state (T « 7000K) at a constant quenching rate of 
2.3 X IQi'^K/s. 

The principle of the thermal conductivity determination 
is illustrated in Fig. 1 . We consider two plates P_ and P+ 
perpendicular to the Ox axis and located at x = — L/4 
and X = +L/4. These plates have a width 26 along Ox 
and their surface is L^. The positions of these plates 
permit to keep the periodic boundary conditions without 
introducing an asymmetry in the system. This has the 
advantage, compared to other studies in which the 
introduction of the thermostatic plates breaks the sym- 
metry, to use a relatively small number of particles. At 
each iteration the particles which are inside P_ and P+ 
are determined and their number is respectively and 
N^. Once these particles are determined a constant en- 
ergy Ae is subtracted from the energy of the particles 
inside P_ and added to the energy of the particles in P+ . 
By imposing the heat transfer in this manner we insure a 
constant heat flux per unit area Jx |^6j which is equal to 
Ae/(2L^At). (the factor 2 comes from the fact that the 
heat flux coming from the hot plate splits equally into 
two parts to reach the cold plate). The energy modifi- 
cation is done by rescaling the velocities of the particles 
inside the plates. Nevertheless to avoid an artificial drift 
of the kinetic energy this has to be done with the total 
momentum of the plates being conserved. For a particle 
i inside P_ or P-|_ the modified velocity is given at each 
iteration by 



VG + a{vi - vg) 



(1) 



where vg is the velocity of the center of mass of the 
ensemble of particles in the plate and 



a = W 1 ± 



Ae 



(2) 



depending on whether the particles are inside P+ or P_ . 
The relative kinetic energy is given by 



rriiVG 



(3) 



Following the standard definition of the transport co- 
efficients 1 16 the thermal conductivity is given by 



Jx 



dT/dx 



(4) 



where dT/dx is the temperature gradient along Ox. This 
formula, known as the Fourier's law of heat flow, is only 



valid when a stable, linear temperature profile is obtained 
in the system. To calculate the gradient we divide the 
simulation box into Ns "slices" along Ox in which the 
temperature is calculated at each iteration. Due to the 
periodic boundary conditions we can concentrate only on 
the Ns/2 slices between x = — L/4 and x ~ L/A and have 
a better determination of the temperature in these slices 
since by symmetry arguments these slices are equivalent 
to the Ns/I slices located outside [—L/4, L/4]. Wc can 
therefore determine the temperature Ti {i — 1, ..Ns/2) of 
each slice at each iteration. By averaging each Ti over a 
large number of iterations to kill the unavoidable large 
temperature fiuctuations (due to the small average num- 
ber of particles in each slice), we are able to determine 
after which simulation time t the averaged profile of T(x) 
can reasonably well be approximated by a straight line. 
After that time we estimate T{x) using a first order least 
square fit of the averaged T^'s, the slope of which will 
give us the temperature gradient. At that point all the 
quantities necessary to calculate k are determined. 
Concerning the "practical details" of the simulation we 
have checked that the results are independent on the 
choice of Ae and for the other quantities we have used 
a compromise between computer time and accuracy of 
the results. Here are the values used in our simulations: 
the width of the plates has been taken equal to 26 = lA 
which means that approximately 30-40 atoms are inside 
the plates at each iteration. The temperature gradient 
has been determined on Ns/2 — 6 slices, each slice con- 
taining approximately 100 particles, k has been deter- 
mined on samples which have been saved all along the 
quenching procedure and therefore have different tem- 
peratures T. To have the same treatment for each sam- 
ple we have fixed Ae to 1% of ksT which appears to be 
a good choice. The temperature gradients obtained this 
way are small enough to insure the validity of Eq. 4. The 
most problematic choice is the simulation time r. Indeed 
in order to reach the steady state one needs long MD 
runs. For us a typical run consists of 50000 MD steps 
(35 ps) directly after the quench during which the aver- 
age temperature is fixed and the heat transfer is switched 
on. Then we perform 450000 supplemental steps (315 ps) 
with only the heat transfer but no other constraints dur- 
ing which the results are collected and averaged. After 
this time the temperature gradient should have converged 
and the value of k should be constant. As we can see in 
figure 2, this can be considered to be qualitatively true 
for the samples above IGK but certainly not for the low 
temperature systems. In fact at low temperature longer 
runs (1 million steps (700 ps)) are necessary and still the 
convergence is not perfect (it is interesting to note that 
though our method converges slowly, it still converges 
faster than the calculation of K{t) given by a steady state 
experiment without a temperature gradient ( p^ , p. 61)). 
It is also worth noticing that the characteristic sigmoidal 
shape of the temperature profile observed at IK is con- 
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sistent to what is expected in the intermediate regime 
where only heat transport over a smaU distance close to 
the plates is effective. In the following, only the results 
above 8K will be reported. 



III. RESULTS 

The results obtained for the thermal conductivity as 
a function of temperature in our model silica glass are 
reproduced in Fig. 3 and compared to experimental data 
cohected between 1 and lOOK and up to lOOOK 
The first observation is that our simulations with the 
BKS potential give the correct order of magnitude over 
the whole temperature range (except at very low tem- 
peratures) with no adjustable parameters apart from the 
"technical parameters" described above and the consti- 
tutive potential parameters. At very high temperatures, 
say above 500K, one observes a more marked saturation 
of k{T) than in the experiments. This might be explained 
by the fact that other contributions than the one de- 
scribed here can occur in the experiments at such high 
temperatures. It is known that the radiative contribu- 
tions (photon transport) in particular increase quickly in 
this temperature range and can become of the order of 
the phonon contributions |l9t . In a large intermediate 
range, 20K to 400K, the agreement between the calcu- 
lated and experimental values is very good. Indeed in 
the simulation also k increases in this temperature range 
unlike what is found in crystalline samples. The major 
discrepancy between the simulation and the experiment 
can be seen between 8 and 20K since we do not find the 
characteristic plateau in the thermal conductivity. In the 
following, we would like to argue that this discrepancy is 
essentially due to finite size effects. 

In our cubic finite simulation box with periodic bound- 
ary conditions, the components of the k wavevectors 
take discrete values of the form k^ = nx^ir/L, where 
is a relative integer (and similarly for the other space di- 
rections), and one cannot find, in principle, propagative 
phonons with a frequency smaller than a lower cut-off Uc 
which can be estimated by 2'kvt / L, where vj- is the trans- 
verse sound velocity. Considering the experimental value 
vt = 3.75 X 10^ cm/s for silica |^ this gives uJc/^-k ~ 1.5 
THz (in practice, when diagonalizing the dynamical ma- 
trix in our low temperature sample, we find, similarly to 
a previous work done on the same system a slightly 
lower first non-zero frequency LOojl-K ~ 1.2 THz, in agree- 
ment with the existence of an excess of modes (maybe 
non-propagative), the so-called Boson peak ||2^, in this 
frequency range . Therefore using the correspondence 
— SkgT which gives the average phonon frequency uj 
of the phonons excited at temperature T, there are cer- 
tainly not enough phonons excited at temperatures be- 
low To ~ 19K in our box to be able to reproduce the 
experimental curve correctly. In Fig. 3, the departure be- 



tween our simulations and experiments is actually seen 
at a temperature of the order of 20K, in good agreement 
with this analysis. 

To try to put this argument on more quantitative 
grounds, let us assume that the thermal conductivity is 



given by the usual formula |24|, 



-Cvi 



(5) 



where C is the heat capacity per unit volume, v and £ the 
velocity and mean free path of the phonons, respectively. 
When applying such a formula to glasses one has to be 
careful because of localization effects. Obviously v and i 
are the characteristics of the "propagative" phonons, i.e. 
those which really contribute to the transport phenom- 
ena. Consequently the heat capacity C to be considered 
should be only due to the contribution of these phonons 
and therefore (according to other authors H,^) should 
exhibit at low temperature the usual Debye behavior (the 
same as in crystals) . If we assume also that the lack of 
phonons in our box, i.e. a wrong value of C, is the es- 
sential cause for the underestimated calculated value of 
K, a very simple and crude way to take care of this is 
to multiply our simulation results by a corrective factor 
Coo/Cb which can be estimated by taking for Coo and Cf, 
the heat capacities calculated in the Debye approxima- 
tion for an infinite system and a finite cubic box of edge 
L, respectively. To calculate this temperature dependent 
factor we have used the standard formulae |24l 



27r2 



Cb 



2 



huj/2kBT 
smh{huj/2kBT) 



hvpk/2kBT 
sm]i{hvpk /2kBT) 



^du (6) 



(7) 



with Lul, = {N/L^)18n^{l/vl + 2/v^) . In the expres- 
sion of Cb the double sum runs over the three polariza- 
tions p — L,Ti,T2 and over the first TV k vectors (quan- 
tized as indicated above) of lowest norm k = \ k \ . For N 
and L we have taken the simulation values N = 648 and 
L = 21.48A and for the sound velocities the experimental 
values vl = 5.9 x 10^ cm/s and vt^ = = 3.75 x lO'^ 
cm/s [20| ]. When correcting our numerical data this way, 
we obtain the open squares represented in Fig. 3 which 
turn out to be in very good agreement with the experi- 
mental results in the plateau region. Of course, our rea- 
soning is very crude since it assumes that finite size cor- 
rections affect only the heat capacity contribution in the 
expression of k (Eq.5) and that the harmonic approxima- 
tion holds for the propagative phonons in that temper- 
ature range, however we think that the agreement with 
the data cannot be fortuitous. It is unfortunate that 
we could not obtain more reliable results at tempera- 
tures lower than 8K (due to the impossibility to reach 
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the permanent regime). Anyway, after correction, these 
results would certainly give larger values for k than the 
experiments since it is known that, at very low temper- 
atures, the propagative phonons start to be scattered on 
the quantum two level systems Q and therefore should 
have a lower mean free path than the one obtained in a 
classical calculation like the one performed here. 



IV. CONCLUSION 

In conclusion we have presented the results of an exten- 
sive classical molecular dynamics simulation aimed to de- 
termine the thermal conductivity in a model silica glass. 
This determination has been done directly inside the MD 
scheme with the use of the standard equations govern- 
ing the macroscopic transport coefficients and no pre- 
conceived model has been assumed. Moreover it turns 
out that this method has considerable advantages (espe- 
cially concerning the length of the simulations) compared 
to the standard methods usually implemented to calcu- 
late the transport coefficients . The calculated values 
of the thermal conductivity are in good agreement with 
the experimental data at high temperature (T > 20K) 
and by including finite size corrections in a simple way 
we are able to reproduce the plateau in the thermal con- 
ductivity around lOK, which has been the topic of several 
interpretations in the literature The agreement be- 
tween the calculated and the experimental values of the 
thermal conductivity is even more striking when taking 
into account the ultra-fast quenching rate used to gener- 
ate our amorphous samples. This shows once more the 
good quality of the BKS potential which permits to re- 
produce the thermal anomalies of vitreous silica with no 
additional parameters. 

Of course, our arguments on the finite size effects should 
be tested in the future by running larger samples. Nev- 
ertheless the simple phenomenological correction is so ef- 
ficient that one can reasonably claim that the initial dis- 
crepancy between the calculated and experimental values 
of the thermal conductivity is indeed due to finite size ef- 
fects and not to a weakness of the method. Therefore we 
believe that this technique is a good way to calculate the 
thermal properties of materials directly inside molecular 
dynamics simulations. 

Most of the numerical calculations have been done on the 
IBM/SP2 computer at CNUSC (Centre National Univer- 
sitaire Sud de Calcul), Montpellier. We would like to 
thank Claire Levelut and Jacques Pelous for very inter- 
esting comments. 
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FIG. 1. Schematic representation of the method used 
to determine the thermal conductivity. More details can be 
found in the text. 
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FIG. 2. Values of the temperature as a function of x in the 
slices located between x = —L/4 and x — L/i for 4 different 
samples and the corresponding least square linear fit: 
(a) T IK; (b) T ^ IIK; (c) T ^ 27K and (d) T ^ 89K. 
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FIG. 3. Log- log plot of the thermal conductivity as a 
function of temperature in silica: 

• : experiment; *: simulations; □: simulations with finite-size 
corrections. 
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